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Abstract: Formulation is given for efficient parabolic equation solution of radiowave 

propagation in inhomogeneous atmosphere and over irregular terrain. Both standard 
and wide angle parabolic equation derivations are presented. Impedance boundary 
conditions are used to characterize the ground. A tropospheric boundary condition 
based on the exact solution of Schrodinger equation in a quarter plane is derived. To 
permit efficient modeling of the irregular boundary, the parabolic equation together 
with the boundary conditions is transformed into a numerically generated curvilinear 
coordinate system. Finally, formulation is presented for a finite difference solution 
using Crank-Nicolson implicit scheme. 
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1. Introduction 



It is well known that multipath fading can significantly effect the link reliability in a 
communications system or target detectability in the case of radar. The path between 
a transmitted and a receiver is often obstructed by natrual or man-made obstacles 
such as hills, buildings, atmospheric layers, trees, rain, fog, etc. Propagation outage 
due to multipath fading depends in a complicated manner on propagation climate, 
terrain features, path length, radio frequency, and fade margin. In the case of atmo- 
spheric multipath fading, interference due to two or more super- refracted rays arriving 
at the receiver via different paths can lead to a complete loss of signal. Other patho- 
logical phenomenon such as obstruction fading (caused by sub-refractive atmospheric 
effects) and ducting (caused by extreme super-refractive effects) are also possible. 
Such phenomenon are more common in warm and tropical climates, particularly near 
shores, where elevated inversions are formed easily due to the large temperature and 
partial pressure differentials. Reflection multipath fading, which is due to interfer- 
ence between the direct and the ground reflected ray depends strongly on the terrain 
geometry and ground constants. Moreover, elevated terrain features could completely 
mask a receiver from a transmitter leading to severe loss of signal (in some cases, it 
is advantageous to site antennas behind hills to provide shielding against undesirable 
interference). It is very important to assess the effects of environment on the link. A 
computer model that can take into account a given refractive index profile, terrain 
elevation data, and varying ground parameters will be very helpful in predicting the 
link performance. 

In this report we present formulation details for an efficient numerical solution of 
wave propagation in an inhomogeneous atmosphere and over irregular terrain using 
parabolic equation. 

Unlike all other previous formulations of the parabolic equation, we will use a modified 
Helmholz equation for propagation in an inhomogeneous atmosphere as suggested by 
Maxwell’s equations (all previous formulations use a Helmholtz equation which is 
only true for fields in a homogeneous medium). Because the parabolic equation is a 
full-wave method, it will include all aspects of wave propagation such as reflection, 
refraction, diffraction, and surface wave propagation. In this respect it is far superior 
to the commonly used ray method. 

Parabolic equation approximation to an elliptic partial differential equation, which 
the true fields satisfy, has proven to be a viable approach for studying propagation 
problems in underwater acoustics. The method is just gaining popularity with the 
electromagnetic community. Although the parabolic equation regards waves as es- 
sentially traveling one-way, it allows a rapid solution of the fields by way of marching 
along the range starting from an initial range. Another advantage of the PE method 
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compared to the ray methods is that it is valid even in the shadow region where the 
simple ray methods completely break down. Furthermore, it appears to be the only 
practical method for predicting propagation over long ranges (greater than 1 km) over 
a wideband from HF (a few MHz) through SHF (a few tens of GHz). The method is, 
however, not without limitations. In its standard form, the accuracy of the method 
is limited to waves traveling essentially within ±10° from horizontal. Furthermore, 
treatment of the boundary conditions on the uneven terrain is difficult. 

The method we propose to use will attempt to remove both of these deficiencies. 
Firstly, we use a Helmholtz-like elliptic equation to describe the fields and arrive 
at a wide-angle parabolic equation subject to certain approximations. To facilitate 
imposition of the boundary conditions on the irregular terrain, the equations will be 
transformed to a body-fitted curvilinear coordinate system. The PE will be solved 
using finite differences on a non-rectangular mesh. This is a major departure from 
previous approaches which will not only make the method more efficient but also 
more accurate. 

In Section 2, we derive the exact equations satisfied by 2D fields in an inhomogeneous 
atmosphere. Impedance boundary conditions are used to characterize the ground. In 
Section 3, we present details on the impedance boundary conditions. Starting from 
the exact equations presented in Section 2 for the fields, we derive, in Section 4, a 
parabolic equation (PE) valid for narrow angle propagation. This case is termed as 
the standard PE. The standard PE is generally valid for propagation angles that 
are within ±10° from horizontal. To accomodate waves traveling at larger angles, we 
present the derivation of a wide angle PE in Section 5. To truncate the computational 
domain we derive boundary conditions on an upper boundary, which are termed as 
the tropospheric boundary conditions. The derivation is based on the solution of 
Schrodinger type parabolic equation in a quarter plane x > 0, y > 0. This is presented 
in Section 6. 

For an efficient numerical implementation, we transform the differential equation and 
boundary conditions to a curvilinear coordinate system. This is presented in Section 
7. Details on the numerical generation of the curvilinear coordinate system are given 
in Section 8. Finally, in Section 9, we present steps leading to a finite difference 
solution of the equations using a Crank-Nicholson implicit scheme. 
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2. Solution of 2D Fields in an Inhomogeneous Medium 



Consider an electric source producing fields in an inhomogeneous region as shown 
in the figure below. Let us assume that both the sources and the medium are two- 
dimensional in nature in that all quantities are independent of the z-coordinate. As 
in the case of a homogeneous medium the fields can be decomposed into a TE Z case 
(vertical polarization) and a TM Z case (horizontal polarization). It is assumed that 
propagation takes place in the x?/-plane. 




From Maxwell’s equations, we have 

Vx£ = -junH (1) 

V x // = jucE + J (2) 



In the case of vertical polarization, the fields could be written in terms of the z- 
component of the magnetic field, H = z// z , ( T E z fields). Substituting into (2) we 
have 

V x H = V x (zH z ) = V// z x z = ju>cE + J => cE = — : - 



In a source- free environment, we have 



cE = 



Vff z x I 



Substituting into (1) we get 



V x 



E = 









= 
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The equation satisfied by the magnetic field is then 






— V , 

ju V e 



) = —jufiH z 



v - Or 1 ) + u2fiHz = ° 



Vertical Polarization 



Once the magnetic field is determined the electric field is given by 



^ _ z x Vi/, 

jut 



Note that H z does not satisfy the Helmholtz equation unless t is constant. 
For horizontal polarization on a similar analysis with E = zE z shows that 

Horizontal Polarization 




( 3 ) 

( 4 ) 

( 5 ) 



If the medium is non-magnetic, y = yo and E z satisfies the Helmholtz equation. We 
will characterize the ground in terms of impedance boundary conditions [3]. 

We may combine the vertical and horizontal cases shown in (3) and (5) into an 
equation of the form 

V • (aVt/>) + flip = 0 (6) 



where 



a = 



P 



M 



TE or Vertical Pol. 



e 

1 

V 

2„2 



— TM or Horizontal Pol. 



Q^n 2 (x, y) > 0 
f H z TE Pol. 



( 7 ) 

( 8 ) 
( 9 ) 



E z 



TM Pol. 



The quantity n(x, y) is a position dependent refractive index of the medium. The 
partial differential equation (PDE) given in (6) is elliptic, for we have 

d ( dip\ d ( dip\ 

& r &) + ^( Q a 7) +/? ' ! ’ _0 
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or 



a 



5 2 0 5 2 0 



50 da 50 5a , 

+ -z 1 - -k~ + /?0 = 0 

51 dx dy dy 



dx 2 5j/ 2 

The discriminant for the above PDE is — 1 < 0 implying elliptic nature of the equation. 
Equation (6) could be expressed as 



1 da 



1 da 



0 



V 2 0 + — — 0 r + - T" 0y H 0 — 0 



a 



dx 



a dy 



a 



For a non-magnetic, loss-less medium, we have 

1 



a = < 



Vertical Pol. 



— , Horizontal Pol. 



CqC 

1 

1*0 



so that for vertical polarization 



1 da 
a dx 



d (\ 



1 de r _ 1 5n 2 

T dx\t T ) t T dx n 2 dx 

2 dn . 

= — (n is the refractive index) 

n dx 



Similarly, 

Letting 



1 da 
a dy 



2 dn 
n dy 






a 2 (x,y) = < 



2 dn 
n dx 

0 

2 dn 
n dy 

0 



Vertical Pol. 
Horizontal Pol. 
Vertical Pol. 

Horizontal Pol. 



equation (10) may be expressed in the form 

V 2 0 + Qj 0x + (Z20y 4- A?Qn 2 0 = 0 



( 10 ) 



( 11 ) 
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3. Impedance Boundary Condition 



Impedance boundary condition relates the tangential components of electric and 
magnetic fields at the interface of two media. If v is a unit normal and s is a unit 
tangent as shown in the figure below, the boundary conditions are given by [3] 




A 

s 






A 

5 



\ 



s = x cos 6 + y sin 6 
v = —x sin 9 + y cos 6 

x — s cos 6 — v sin 6 
y = s cos 6 + 0 cos 0 

I (12) 

where A s = Z s /r/ 0 is the surface impedance normalized to the free space value T] 0 . 
The equation may also be expressed as 



(z> • E)0 — E = —T] 0 A s i> x H => v x E = tj qA s {/ x (D x //) 

= tj q A s [(i> ■ H)v - H 



(13) 



The surface impedance is determined from the intrinsic impedance of the medium 
by considering plane wave reflections from the interface. The complex propagation 
constants, 7 1} 72, and the intrinsic impedances rji and 7? 2 i n terms of the media 
constants are indicated in the figure. 




7i = 



72 



+ jue 0 £rl) 

-koUrl Url ~j— ) 
^o/J-r2^-Tc2 




According to Snell’s law, we have 7 j sin#, = 7 2 sin# ( 
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The plane wave reflection coefficients for the vertical and horizontal polarizations, Ry 
and Rh are [7] 

t/ 2 sec 6 t — Tj] sec 6 X 



Rh 

Ry 



r\ 2 sec 9 t + r/j sec0 2 

tj 2 cos 6 t — Tj\ cos 6 X 
t) 2 cos 6 t — 7/1 cos 6 X 



Surface Impedance = 7/ 2 sec0 t 
Z\ — 7/2 cos 0 t 



A" = 



7/2 sec 0 t rj 2 
Vi 



1 - lysinfl, 



- 1/2 



//7 r 2^rcl 


1 Hr\£rc\ 2 / 

1 COS XPi 


1 / 


y 1^t\ €tc2 


^t2^tc2 



1 - 1/2 



and 



I^t2^tc\ 


(l T \ t TC \ 2 / 

1 COS Xpi 


J 


y 1^t\^tc2 


^t2^tc2 



Ar = 

For the special case of /j,\ = // 2 = Ho- — 0, 



1/2 



Af = 



a; = 



For normal incidence ijj t = 90° 




1 - 1/2 



1 1/2 



(14) 

(15) 



1 



Af = A r = J — 

V Cr2 - ]°r1 

For the 2-D case, impedance boundary conditions for vertical polarization can be 
simplified as 

0 x E = 770 A]f [(*> • H)v -H] = -rjoA^'H 

Taking a dot product on both sides with z 
z ■ (0 x E) = -Vo^'H z 

Since z x 0 = — s, we have 

-s • E = -Vo^H z 



Substituting from equation (4) for E , we get 

z x V// 



JU>C 



= -VoA V S H Z 
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v ■ VH Z = jueqoA^ H z = jk Q e T A v s H z 
i.e., the above can be put conveniently in the form 



dH z 

dv 



— jk 0 e r A; H z = 0 



IBC Vertical Pol. 



A similar analysis for horizontal polarization yields 

IBC Horizontal Pol. 



dE z 1 „ n 



This may also be obtained by resorting to duality. 

For a perfectly conducting material, we have A^' v = 0 and 




0 Vertical Pol. 

0 Horizontal Pol. 



( 16 ) 



(17) 
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4. Standard PE Derivation 



We will make some approximations and cast (11) in the form of a parabolic equa- 
tion which permits a rapid numerical solution. 



O i^QX 



Let rp{x,y) = — r^u(x,y) 



\/x 



Then 



/ u \ e~ ]k ° x e~ ]K ° 

e ~jk ox e -jk 0 x 

Xjly = Uy ~=~ , Uy X = U T y ~ (U Xy ~ jk Q Uy) ‘ 



o -jk 0 x 



x — * oc 



\/? 



y/x 



V’yy u y 



?~]k 0 x 



? / 

U'.r.r 



(u rr - 2jk 0 u x - klu)- 



o~jk 0 x 



y/x 



VV yfi ’ 

Substituting into (17) we get 

u xr + u yy - 2jk 0 u x + - 2jk 0 aiu + a 2 u y + (n 2 - l)^u = 0 

or 

U xx + Uyy + (a) - 2 jk 0 )u x + a 2 u y -f ^n 2 - 1 - 2 'jj^j klu = 0 

If now we impose the approximation that 

l^xxl ^ l u x| ) 



we obtain 



-1 



Ur- = 



(aj - 2 jk 0 ) 



Uyy + a 2 Uy + ^7i 2 ~ \ ~2j ^ ^ U 



or 




(18) 



This is the exact form of narrow angle PE approximation. We would also like to 
express the impedance boundary condition in terms of the l u ’ functions. 




A 

s 



x = s cos 6 — v sin 0 

dx , „ „ 

x„ = — = v • Vi — v • x 
du 

— — sin 6 

v = — x sin 0 -f y cos 0 
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For vertical polarization we have 



-^--jk 0 t T A v s H z = 0 =► -^--jk 0 n 2 A v s H z = 0 



e~ jk 01 



H z {x,y) = 


u(x,y) 

v X 






dH z 


— jk 0 x u u 


ux v ^ 


e~ jk ° x 


dv 


2x 3 /2 ) 








e -jkox 






(u„ — jk 0 x„u) 


VT ’ 


x — * 



■ ■ u u — jk 0 (n 2 Al' + x u )u = 0 
Similarly for horizontal polarization, we have 



IBC Vert. Pol. 



u v - jk 0 




IBC Horz. Pol. 



We combine the two by defining 



' — jk 0 (n 2 A] — sin#) Vert. 

~j k 0 ( ~ Sin 0 J HorZ ' 



and writing as 



+ Cju = 0 



(19) 



The parabolic equation given in (18) is valid for propagating angles close to horizon- 
tal (±10° in practice) [1]. To accomodate waves at higher angles we would need a 
wide angle parabolic equation whose derivation is accomplished through a pseudo- 
differential operator formalism [1]. 
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5. Wide Angle PE Derivation 



Let us assume that the media constants are independent of horizontal range so that 
a(x,y) = a(y), n(x^y) = n(y). We then have from (17) 



Let 



u T 



+ u 



yy 



-j2k 0 u x + — + ( n 2 (y ) - 1 )k^u = 0 

a ay 



P 



d Id 2 1 da d 

dx ’ \ kl dy 2 kfia dy dy 



( Pseudo differential 
operator in y 



Using this notation, the above PDE may be expressed as 

P 2 - 2jk 0 P + (Q 2 -l)k 2 o }u=0 

which we may factorize as 



{P ~jk 0 — jk 0 Q)(P -jk 0 + jk 0 Q)u = 0 

To see this we expand the operator on the left hand side to get 

P 2 - jk 0 P - jk 0 QP - jk 0 P - kl - klQ 
+ jkoPQ + k 2 Q 2 + k 2 Q Q 



( 20 ) 



Since PQ = QP, we have the desired result. The first operator in (20) denotes 
an incoming wave (w'.r.t. x) and the second an outgoing wave. We retain only the 
outgoing wave to obtain 

Pu = -jk 0 (Q - l)u (21) 

The square root operator Q is global in nature and we would like to make some 
approximations to derive a local operator from it. (It is global because when expanded 
in terms of series, it will contain terms of all derivatives). Now 



Q = 



11 da d Id 2 
kl a dy dy + kl dy 2 



1/2 



Let us rewrite Q as 



Q = 



1/2 



i + W(y) — i] h — 



1 da 



+ 



d 2 



<<1 normally 



a d(k 0 y)d(k 0 y ) d(k 0 y ) 2 



small 



12 



which may be expressed as 

Q = vTTu 

where 



V 



n 2 — 1 + 



1 da d 
k 2 0 a dy dy + 



H dy 2 



a small operator! 



Treating V as an algebraic factor < 1, we may derive the following rational approxi- 
mation (pade( 1 , 1)) 



Q = vT+U 



so that 



nir 

1+ i 1 ' 

4 + 3V 
4 + K. 

2V 



(Claerbout) 



= 1 + 



2V 

4 + V 






Substituting into (21) we arrive at 

Pu = jk 0 (Q - l)it 



—2jk 0 V 

P u = — — u 

4 + K 



or 



Using the fact that a 2 = — 

a 



(4 + V)Pu = —2jk 0 Vu 
d a 

— , we get 

dy 



/ 2 OM 2 ^ d 2 

(n + 3)t 0 + °-*f y + Qji 



u x = -2 jko 



(n 2 — l)^o T <2 2 '^ — h 



d 2 



dy dy 2 



( 22 ) 



The above equation is a wide-angle parabolic equation valid for propagation up to 
±20° [1]. Other approximations could be obtained by considering higher order pade 
approximants. 
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6. Boundary Condition on the Upper Boundary 



To truncate the computational domain, we consider a point high enough where the 
atmosphere is homogeneous with n = 1 . The governing equation in the homogeneous 
region becomes 



u T — — - 






2*o yy 

Let us derive boundary conditions on a horizontal interface y = y 0 . For the sake 
of simplicity we will derive boundary conditions of the mixed type on the interface 
y = 0 (instead of y = y 0 ) given the initial data on the line x = 0 and boundary data 
on the line y = 0. Our derivation is based on the use of Fourier sine transforms as 
suggested in [2]. Although the basic philosophy of our approach coincides with that 
in [4], some of the details and the final results are slightly different from the latter. 



Consider the parabolic equation u yy — 2 jku T = 0 in a homogeneous region x > 0, 
y > 0, where A; is a complex constant, 




subject to the initial condition u(0, y) = /(y), 0 < y < oo, and the boundary condition 

u(x,0) = g(x), 0 < x < oo. We assume that u(x,oo) — * 0 and u y (x, oo) > 0. The 

equation 

u yy = 2jku x (23) 

is of Schrodinger’s type. We will treat the lossless case having a real value of k as 
the limiting case of the lossy problem having k = k 0 — je, t > 0. We will solve the 
problem using Fourier sine integral. 



Let 



l~2 /*°° 

U s {x, A) = \ - / u(x,y)sm\ydy 
V 7T Jo 



(24) 
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Then using integration by parts, we see that 




7 r Jo 



u yy{x, y) sin Xy dy = \ u y (x,y) sin \y 



7T 



y=0 



—Xu(x, y) cos X y 



y=o 



roc 

— A 2 J u(x, y) sin Ay dy 



Because u v (x,oo) — ► 0 and u(x,oo) — * 0, we have 
I ~2 r°° 



— u yy (x, y) sin Au dy = \ — Au(x,0) — X 2 U s (x, A) 



77 JO 



77 



= \l -\g(x) - \ 2 U s {x,\) 



(25) 



Multiplying both sides of (23) with \J2/77 sin Ay, integrating over y = 0 to oo and 
making use of (25), we have 



or 



^Xg(x) - X 2 U s (x, A) = 2jk-^-U s (x,X) 



dx u s {x,x) + 2jk u *( x ’X) - 2jk ][l 9 ( x ) 



We may rewrite the above equation as 



d 









(26) 



2 jk V 7 r" 

Replacing the dummy variable x in (26) with r and integrating both sides over r = 0 
to x, we arrive at 



U a (T, X)e^ 2/2]k) ' 



X x 



r=0 2 J k V * Jt =° 




g{T)e ix2/2]k)T dT 



or 



But 



U s (x,X) = U s (0,X)e-l x2/2 ^ x + ^ ' 2 ^ 



\l^ I" g(T)e {x2/23k){T ~ x) d7 
V 7T Jr= 0 



f/ s (0,A) = J- f u(0, y) sin Xy dy 

V 7T JO 

= <F /"/(</) sin Xydy = F,(\) 

V 77 Jo 
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u.(x, A) = F ,( + -2-J- r g{T)el*n k ^-*UT (27) 

2] K V 7T 7t=0 

Finally, taking the inverse sine transform on (27) we get 

u{x,y) = \/- / sin Ay t/ 5 (x, A)dA 

V 7T Jo 

= \h- f°° F s {\)e- {x2/2]k)z sm\ydX 

V 7T Ja=0 

+ --TT /°° A sin Ay /* y(T)e (A2/2jfc)(T - x) dr dA 
7r2jkJx=o Jt = o 



7 r OC rOO 

= — / / /(r) sin At dre~^ x sin Ay dA 

7T 7a=o A=o 

+ --^- r y(r) r \ sin Xye {x3,2ik){r - x) dX dr 
7r2jkJ T =o J \= o 

1 /* oc /*oo 

= — / /(r) / [cos(r — ?/)A — cos(r + y)A] e”^ A ^ 2jk ^ T d\ dr 

7 r Jr=0 Ja=o 

+ - f g(r)- T (-f) [°° cos Xye^/^-^dXdT 
7 r Jr=0 J A: \ oy ) J A=0 

1 /*°° r°o 3 

= — /(r) / [cos(y — r)A — cos(y + r)A] e”^ A Fi k ) x d\ d T 

7T Jr=0 dA=0 

-^1 JL 9{T) Ty {/„” C0S!,A e ' (i,/2 ’‘ ,( '" ,<iA } * 



Defining 



A'(x,y; x 0 ,yo) = [ cos (y - y 0 )Ae ( aJ /2jD(* for x > X()j 

./o 

we write the expression for u(x,y) as 

1 f 00 

u(x,y) = - f(r) [A'(x,y; 0,r) - A'(x,y; 0, -r)]dr 

7 T Jt =0 

1 /* X 5 

— r / 2( T )'a- A '( ;,: >2/; T > 0 ) rfr 

7ryK i/ T— o Oy 

We now evaluate the integral for A'. Consider 

/(a) = / cosQAe- (A 2 / 2 jfc)(x - Io) dA, x > x 0 

7a=o 



( 28 ) 
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= r cosaAe (A2/2| * |2)j<: * (x " Xo) dA 
J A=0 

= r cosa\e-( x2{x - xo)/2m <- jko) d\ 
J A=0 



Because of the exponential decay, we may differentiate under the integral sign to 
obtain 

J-I(a) = - [°° AsinaAe- (A2(l - xo)/2|,:|2)(£ - jio ^A 
da Jx= o 

e -(A 2 (i-r 0 )/2|fc| 2 )(£-y<:o)l 



= J s\ 

Jx = 0 



sin aX 



d_ 

d\ 



(x - x 0 ){e -jk 0 )/\k\ 2 



dX 



sin Q Ae- A2 ( I - IO H t - 2 fc o)/ 2 |<:| 2 



a|Af 



(x - x 0 ){e - jk 0 )/\k\ 2 

I™ cos QAe- (A2(x - Io)(£ - jJco))/2| * |2 (fA 
Jo 

jo '\k\ 2 



=0 {x - x 0 )(e- jk 0 ) 



(x — x 0 )k 



-■L 



°° cos aye- x ^ x - X0 ^ 2jh d\ = - - J — - 1(a) 



(x - x 0 y 



dl(a ) jak T . . 

4r + <^) /(Q) - 0 



or 



= 0 => /(a) = I(a = O)e- J “ 2t/(2(I - Xo) ) 



— \l(a)e J<y2k ^ 2 ( x - xo ^ 
da I* v 

Now I(a = 0) can be written as 

I(a = 0 ) = 1°° e - (x2(x - Xo),2 ^ )(t - 3ko) d\ 

Jo 

Let us view this integral in the complex A-plane. 




(29) 



I 



For the integral to converge for ( x — x 0 ) > 0, Re[\ 2 (t — jk 0 )] > 0, i.e., 

i?e[A 2 (— >A:“)] > 0 Im(\ 2 k’) > 0 



or 



0 < Arg (A 2 /:*) < 7r 



or 



or 



Now 



0 < 2 Arg (A) — Arg ( k ) < n 



\ Arg (k) < Arg (A) < j + ^ Arg (k) 



Arg (k) = — tan 



-l 



ko 



ko ko 



< 1 



6 A /XA ^ 6 

•• ~W 0 < A,gW < 2 - W 0 

I(a = 0)= / e~ x ' 1(x - xo ) ( - ]k ')l 2 \ k? d\ 

Jc 

where C lies in the region of convergence. 

In particular, let us choose a line from 0 to oc along the line C w defined by 



(30) 




On this path, 



A / \ \ 1. / 1 \ ^ 

Ar6(A) = j+-Arg(i-) = 7 - — 






TTjk 



2(x - Xq) 
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/(a) = ./ * jk 

V 2(x — x 0 ) 






( 32 ) 



Substituting this into (28), the field u(x,y ) for the initial-boundary value problem is 
given by 

u(x,y) = - /~ /(r)\/^{e- J ^- T)2/(2r) -e-^ +T ) 2/ ( 2l )}dr 

7T Jr=0 V 2x ^ J 



■ _L f g( T )—l.j *> k I e -^/W«-r))l rfT 

njk Jt=o dy \\j — t) J 



It is easy to see that 

f = 

d 2 I< 



dy 2 



= /°° -A 2 cos[(y - y 0 )A]e- A2(r - T)/(2j,c) dA 

Jo 



and, so for x ^ x 0 



<PK_ _ dK_ dK_ 1 d 2 I< 

dy 2 ^ dx dx 2 jk dy 2 



which is the same equation satisfied by u. Now 

u(x,y) = \fp- r /(r) f e -^- T ) 2 /( 2r ) - e - jk(y+T)i/(2x) 
V Z7 TX Jr- 0 

1 r x d 

-T / 9(r)^-I<{x,y, T,0)dr 

7rjkJr=0 uy 



\ liH 





1 jk j 


l 

0 

+ 

1 


V 2nx Jo 



1 TJ 



[~ f{T )^ e -ikr*nu) dr 
Jo X 

rkJL 3{T) h K{x,v: T,0) 

jk r f(T)2^-e~ lkT7/{2x) dT 

Jo OT 



dr 



y— >0+ 



27TX 



r,o) 



njk 



dr 



y— +0+ 



(33) 



(34) 
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in view of (34). Furthermore, we note that 

(x, y, xo, yo) ~dx ^ 1 *^o 2/o ) 

Using this and evaluating the first integral by parts, we get 

d 



dy 



u(x,y) 



y — >0+ 



= -/^{/( r ) e " jfcT2/(2r) - I" fr(T)e- ]kT * /{2X) dT 

2 r* d I 

+ “ / flr(r) — A'(x,y; r,0) dr 

7r Jt=o Ot , 



ly— 0+ 



= {-/("> - J o °° MtK’^t} 



2 

+ - 
7T 



y(r)A:(x,y;r,0) 


T _ 0 “ / 5 r r(r)/\(x, 0 ; r,0)dr 




v-o+ A =° 



= +J—m + 

V 7TX 



I2jk 



7TX 



rOO „ 

/ / T (r )e~ 3kT ^ 2x ^di 
Jo 



-b m ' 



'njk _ 1 2 jk r* y T (r) 

V 7T Jt = 0 \fx~—~T 



2x 



v/3 



dr 



\2jk 



7TX 



im - 9 ( 0 )] 



Wi \ r 

7 r I 7o 



)rfr - f 

Vo v X — T 






\/x do \/3 

From the compatibility conditions on the initial and boundary values, we have 



Therefore 

du 



dy 



(x,y) 



lim u(x, 0) = y(0) = limu(0,y) = /( 0) 

x — ►O y— *•0 



4 = r fr(r)e~ ]kTi/{2x) dT- r -pLLdr 

\JX J t= 0 7r=0 — r 



y— -0+ 




9r( T ) 



( 35 ) 



It is easy to note from (35) that for k = k 0 — je, 



du 



= 0 . 



y— 0+ 
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No matter how small the c is, we will use the same equality for the limiting case 
of e — > 0. We will evaluate the integrals above approximately by replacing the 
derivatives with differences. 



y t 



Ay 



I 



u, 



IL 



U. 



U 



Pvi'rJlS 



\ 



u u 3 



X 



Consider initial data on the line x = 0. Let us assume that this initial data is known 
on a uniform grid y m — mAy, m — 0, 1, - - We approximate the derivative in the 
interval {y m -\,y m ) by the forward difference formula 



^ y ( V ) ^ ^ '7i ^ m — 1 

dy Ay 



y e (j/m — 1 1 Vm ) 



where u m - u(0,y m ). 



Then 



Let 
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We then have 



\A/(’ rr )s' m 

^/l"„.F ,(tT ' ,/2Irfr - vf / 



= \/l F M 



\Jkf (^x)y m — 1 

\A7(^)y m 
y/h/{7rx)y m — 1 



where F(/i) = complex Fresnel integral [8] 
1 fC ° ' ' ‘ c -ifcrV(to) dT = 



1 /*°° 

-f/ A(r) 

^/X J T=0 



E / U m U rn — \ \ 111 

A Ay Jv* 




J/n 




2 /m — 1 



(36) 



Consider now the boundary data on the line ?/ = 0. Assume that we have a non- 
uniform grid 0,xi,X2, • • * ,£a/ = x. On the interval (x m _i,x m ), we approximate the 
derivative as 

9t(t) = 

'Em — 1 

where u m = n(x m ,0). At the origin we have u° = u 0 . 



u m _ u m - 1 



We evaluate the second integral in (35) as 

9t(t) M " m - 1 



/*4±l* = f; p 

Vo VX *7” m=l *^m — 1 J x m—\ y *E T 



= n 



M U m - u m_1 



— i Z-m - 1 



^m-l 



- 2E 



M u m - u m_1 



— i 'Em — 1 



m — 1 



A/ 



(—2\/x — T j 
^ \Ar x rn ^ 



= 2£ 



U m - U m_1 



— j ^m- 1 “1” \/Fw 



(37) 



Using (36) and (37) in (35) we get 






y=0+ 



[ 2 jk fW u m u m _i f r-, f I ^ 

v~lv*s — — v F (te 



2 /m I F 



7TXj\,/ 



" 2 /m — 1 



M 



-2E 



u m _ u m-l 



— I \/lM Fn" ~h A/ — 1 
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Extracting out the m — M term from the latter and writing 



we have 



!<*■»> 



■+ 0 + 



du™_ 

dy 



du^_ + 8jk v™_ 
dy V 7T y/x M - 




1 



F 




-F 




u m — u m 1 



7T m = 1 \/ X M '^rn "P \Z-^M ^m— T 



(38) 



The above equation is the discrete version of a continuous boundary condition of the 



form 


du , N 


(39) 




— + r(x)u = s(x) 

dy 


where 


and 


/ \ 8 J k 1 

r x = \ 

V 7T y/x — X M - 1 


(40) 



'M = 

m — 1 ^ 







it 

2/m — 1 



M - 1 



-\Pr£ 



u m — u m 1 



7T m _j \/^ “1” \/~3~ X m — ] 



and 



F(x) = [ X e- j ^ /2 ^dT 

Jo 

is the complex Fresnel integral [8]. 



(41) 



(42) 



For efficient implementation, the PDE and the various boundary conditions will be 
transformed into a curvilinear coordinate system generated by setting the lower ir- 
regular boundary as a constant curve curvilinear coordinate. 
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7. Transformations to a Curvilinear Coordinate System 



Consider the narrow angle parabolic equation with range dependent refractive index 
(a] ^ 0) given in ( 18 ) 



1 



u T = 



d d 2 

a ' + a ^ + W 



U 



(2jk 0 — O] ) 

together with the tropospheric boundary condition 

*5?/ 1 )^(^mi 2/o) ^ (*^m ) 2/o) 



PDE 



and the impedance boundary condition on the irregular boundary 



*7 7 “ 



!i„ + CjU = 0 




<JL 



^(jwvvC^a y 



We will transform these to a curvilinear coordinate system (£,77) 
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We assume that we have a transformation of the following form 

* = *(0i y = y(Z, r i) 

The various metrics needed in the transformed equations are [5] 

9n = x\ + y\ , 912 = x^x„ + y^r, . 

We then have with yjg = x^y v — x v y^ = x^y v ^ 0 

u x = - y ( u v ) = -^—( Ui y v - Uri yz) 

y9 x (Vr] 

U y = -^(-XrjUt + X f U„) = — 

y/9 y*) 



l d 



u yy ~ 





1 


Vr)T) 


u 


„ , 2 


U T)T) U T) 


y 2 v 


[ Vv J 



Substituting into the PDE we get 

1 < ' \ 1 

(u ( y v - u n y ( ) = j— 



xtyr, 



(2 jk 0 - g j ) 



» . u v yjrn , \ 

Ijlt + 02 T-Ur, H — > 

yn y n yl J 



or 



Xca, 

u ( = Z U + 



(2 jk 0 — ax) 



a 2 y 7)?7 



x t 



+ 



x m 



Jv Vv ) ( 2 i k o-«x) Vn 



Xe 

U y "h /r> • ; \ 2 Ur,TI 

(2jk 0 - ax)y2 



Letting 



&i = 



X^Oj 



(2 j ko — ax ) 



62 = — 



Vn 



bz — 



a 2 - 



x t 



Vw 7 



1 



y* / (2j^ 0 - «i) 



+ y* 



(2j£ 0 - ai)yjj 



we express the PDE as 



u ^ — b\ 11 -p 62^77 + 63 u 



T)T) 



(43) 



The normal derivative onai/= constant line is 



, , % 9 11 yi2 

U^[T] = COnst.J = ./ Uj, Ue 

V 9 \/99n 



\M + yl 



X^rj-y^Xr, ” rJxfTyi{x ( y v -y ( Xr,) 



X(Xr, + y^r, 



H 



(44) 



25 



Variation of an arbitrary vector r along the r\ = const, coordinate is 

r K = x ( x + y ( y = 

The unit vector along the tangent on rj = const, is then 



s = 



a c 



H x ViV 



l Q d yjx\ + y\ yjx J+ 



yi 



Defining 



cos 6 — 



sin 6 = 



H 



\R+ 



y( 



Vi 



+ y\ 



we express (44) as 






U v — 



(x v cos 0 + y v sin 9)u ^ 



y n cos 0 - x v sin 0 ^/xl + yfiyr, cos 6 — x v sin 0 ) 



With these substitutions, the boundary condition at the bottom boundary h 
0 gets transformed to 

(x v cos 0 + y n sin 6) 

u„ , — u £ + ( y v cos v — x v sin v )C\ u = U on rj — (J 



\J x l + 



vi 



For x v = 0 and using \J^\ + y\ = x r J cos 0 , we get 



V 

u v sin 6 cos 0u£ + Ciy^ cos Ou = 0 @ tj = 0 

x i 



Finally at the top boundary, we have 

u r,/y n + ru = s 



or 



A ) "h c(£ m , A Af) — yr} S {£mi N) 



v + C\ U — 

(45) 



(46) 



(47) 
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8. Generation of Curvilinear Coordinate System 



Consider a piece-wise linear ground profile and a horizontal upper boundary 




Tj-N 



... ■■ , ,u . 

~£y\A^i^JL 


/ * 

A' 


f— * * f 

B' 


^ f ' t 

c' 


D 

^ 














~il ,k 

AY) 

V 


ti 

X 

? —r* 


> , 


* — -*■ » 
-7 7—^ 


/ / 4 



> y y—i y / k , y Hi 

A 6 C D 1 



Non- Rectangular, 
Uniform Mesh 



Rectangular, 
Uniform Mesh 



Physical Domain 



Computational Domain 



+JL c — « 



*t.y. 






.y. 




N 




We generate the (x,y) coordinates of an interior point by using a linear interpolation. 
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Letting Ax, = (x , +1 - x,), Ay, = (y ,- +1 - y,), and A£, = ( 6 +i - £,), we write 



x = x * + ~ 6 ) 



y 



-0-S) 



, A*/,-, > 



r\ 

+ N y ° 



(.<(< fi+i 
0 < y < A 



At any interior point, £i < £ < £, + i> the various metrics are evaluated as 
Ax, 



x i = 



Vi = 1 



A6 ’ 



x v = 0 



nJ 



y \ Ay, 



A£, ’ A r 



Ay, f e e ^ 

y° - y. - - 6) 



( 48 ) 



(49) 






At the boundary points, we use the central difference formulas to arrive at 



*f(f = C) = 


x , +2 - x, Ax, + 1 +Ax,- 


(50) 


6+2 — 6 A^, + i + A£,- 


<CS 

II 

II 


L V\ Ay,+i + Ay, 
V n) A 6 + , +A 6 


(51) 



Note that the analytical expressions yield a discontinuous value for these derivatives 
at the boundary points. We use the analytical expressions only to generate grid 
points and use the central difference formulas to arrive at the derivatives w.r.t. £. 
In other words, once the grid points are generated on the lines AA\ BB ', . . . , etc., 
we assume that the space is smoothly connected through the grid points. In the 
numerical implementation using Crank-Nicolson implicit scheme [ 6 ], the metrics are 



needed at the midpoint w.r.t. £, i.e., at £ = £, + A£,/2 and the interior point formulas 
are applicable. For a uniform mesh in the computational domain, A £, = 1, y = y, 
q = 0,1,2, •••,7V 


x t = 


Ax, 


(52) 


Vi = 


(i N ) Ay ' =>■ y«(y + 1 ) y<(?) = N 


(53) 


Vr, = 


NV° 2 ) 


(54) 


y = 


f, q\ (y%+ \+yi\ , q 
(‘ n) { 2 ) + N y ° 




— 


2 + •?:/„ = y<> + (<? A')s/, 


(55) 


so that y(q + 1) - y(q) 


= Hry 
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9. Numerical Implementation by Crank-Nicolson Scheme 

Consider the narrow angle PE together with the boundary condition 

u , £ = biU + b 2 U n + 63 Urjr, (PDE) 

u„ — ( — sin 6 cos 0 ) u« + Ciy n cos Ou = 0 at 77 = 0 (BC2) 

\ x ( J 

Uj, + ry v u = y v s at 77 = N (BCi) 

with 

u,( 0,N)=0 

We would like to implement the above using a Crank-Nicolson implicit scheme [6] 




h 

v 



We use the notation u v q = u(£ p ,r] q ) = u(x p ,y q ) 

The various derivatives assuming = Ay = 1 are 

u ti£p-l/2 lVq) = U q U q 

u n (^P ~ 2 '^} = *7?) ^ 2 *7?) Uj ?(£p> Pi)) 

= \ [<+i - <-l + <n - <-i] 

U V = \ K+l 1 + U 9+l “ ( U £-l + <-l)] 



or 
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1 



u nn = \ K+i “ 2 < + u g-i + - 2 < 1 + <-!] 

Substituting into the PDE, we get 

i / 1 \ u? + u p n ~ l 

= '-'(p-j-?) — 5 — 

, b2 { p ~bi 



(<+l + U q + 1 - <-l - <-l) 



1 



« ^ ( U g+1 2u q + U q-1 + U q + \ ^ U q 1 + U j-l) 



Rearranging the terms we get 
1 (b 



v 



1 



2(y + ‘3K +1 -|i + ^K + -U3-f 



1 f bo 



V 



+sl? + M*G! + |i-fe + Tl«r , + s(*»-T <:l = o 



2 \ 2 



Now let 



Q 



-H 



1 ^ P~ 1/2 

b 3 +~y 



» = ( 6 -l 



9 

P-1/2 



7 = iff.3-^' 



P-1/2 



We then have for p = 1 , 2, • • •, q = 1, 2. • • •, N — 1 



mi ^ 



Ui - 0 + 0K + 7<-i = -ortij;! - (i - /?)<-“ - 7<:{ 



,p-i 



,p— * 



(56) 



However, we will extend the applicability of this equation over q = 0, 1,--/V to 
accomodate the derivative boundary conditions on the lower and upper boundaries. 
For q = 0, we have 

ctu p -(1 + /?)ug + 7 u p _ 1 = -ou ? _1 - (1 - /?)ug _1 - 7 u ^ 1 (57) 

for q = fV, we have 

~ (1 + 0)u P n + 7*4-1 = _Q; *4+i - (1 ~ ^) u n _1 “ 7 u n"-i (58) 
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From the BC2 at the lower boundary we have 



M f> 

"f f" 






7 7 * X P y 7 - ^ - 0 

1 1 v 

1)--' 



I 

I 



,P“1 0 ,P“ 1 „,P „,P ■ 



1 / \ P-1/2 

-till , ti;-wii ( y v . ^ ^ / p P -i\ 

2 b 2 sin 0 cos 0 1 [u% - ttg J 



+ (c,,,c OS *r«fc£ 



= 0 



Multiplying this with 4(7)^ and adding to ( 57 ) 
(a + 7 )u p l - ( 1 + /? - 27 



* y 

Ci y n cos 0 — 2— sin 0 cos 0 
x ( 



U n 



= -(a + 7)1*1 1 - I 1 - /3 + 27 



Letting 



a' = a + 7 



V 

C\y n cos 0 + 2— sin 0 cos 0 
x t 



u: 



p - 1 



P' = P - 27 1/,, cos 0 ( Ci - 2sm ^ j 



/?" = P - 2^y n cos dc,+ 



2 sin 



we may write the above equation as 

a'u P i - (1 + P')u p 0 = -a'u P i~ x - (1 - P")u p Q ~ x 
From the BCi at the top boundary we have 

* * M+ | 



X L. 



f ' * ' * <- N 



M 



N-l 



( 59 ) 
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u 



P-1 

N + 1 



— u 



P-1 

N - 1 



+ 



TV-fl - u 



TV -1 



+ 



r p -f r p 1 Um -f u p 1 



-</»>- 



/V S P + S 

= 2/r,— 



P-1 



Multiplying this with — 4(a)^ and adding to (58), we get 

- (l -f /? + [r p -f r p_1 ]i/^a) + (7 + 4 a)u^_! 

= - (l - /? - (r p + r p_I ]i/ r ,a) u^" 1 - (7 + 4 a)u’ 7 _i - 2 y v (s p + s p_1 ) 

Letting 

A = /? + y^a[r p + r p_1 ] 

7 # = 7 + 4a 

we write the above equation as 

(1 + A)uJ, - -/’i' F v-i = (1 - + I'uS, - .', + 2j„(r f ’ + r p_l ) 



(60) 



We augment the equations given in (56) for q = 1,2,..., — 1 with (59) and (60) for 

q = 0 and TV, respectively to define for q = 0, 1,2, ... ,N. The system of equations so 
defined can be expressed as a matrix equation of the form 



X X X 
X X X 



X X 



- 


u p 0 


- 




„ p 






u l 






p 




- 


. Tv J 


- 



X X 
X X X 
X X X 



X X 





„ p-i i 








p - 1 




< 




„ p - 1 




. U N . 



(61) 



where X denotes a non-zero entry. The tridiagonal matrix on the left hand side of 
(61) can be inverted efficiently to yield a solution on line £ = £ p in terms of the field 
values on the line £ = £ p -i- Equation (61) can be used to march forward in range 
starting from initial data specified on £ = 0. 
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